% NUM_REP = 10;       % The number of Repetitions
% ID_SLICE = 8;       % The ID of slice
% 
% FileName_M0 = '';   % The name of the 'M0' file
% FileName_ASL = '';  % The name of the 'ASL' Images

%% Read all asl data into memory
parameters = divasl_set_parameter('default');
[f_list_com0,f_list_coasl,f_list_mri] = divasl_get_filelist(parameters);

% obj_divasl = cDIVASL(f_list_com0,f_list_coasl,f_list_mri);
obj_divasl = cDIVASL(f_list_mri);



%% Read the template images
% templates should be resliced before using
%  job_reslice_file(true,'',false,'c2rT13D_asl_file.nii','cpwall_nrtasl_file_gm.nii');

f_template_graymatter = 'myProject/template/rc1rT13D_asl_file_on_cpwall_nrtasl_file_gm_nii.nii';
f_template_whitematter = 'myProject/template/rc2rT13D_asl_file_on_cpwall_nrtasl_file_gm_nii.nii';
f_template_shell = 'myProject/template/rc3rT13D_asl_file_on_cpwall_nrtasl_file_gm_nii.nii';

obj_divasl = obj_divasl.setTemplate(...
    'gray',f_template_graymatter,...
    'white',f_template_whitematter,...
    'shell',f_template_shell);   
                                        

%% Get only asl at gray matters
% 
for i = 1:11
    obj_divasl.dDivasl(i) = obj_divasl.dDivasl(i).calCBFImg();
end
% If you want to calculate the cbf in white matters, just use the following
% code
% obj_divasl.dDivasl(i).calCBFImg('buxton','mask','white');


obj_divasl.dDivasl(6).dCBF.showMRI('Mono TI')

%% Compute the CBF by diversity ASL, using Least Square Method(dCBF_LS), 
%   Compressive Matched Filter(dCBF_CMF)
obj_divasl = obj_divasl.callLS('Gray','all',1:2:11);
obj_divasl.dLS.showResult('cbf',obj_divasl);
obj_divasl.dLS.showResult('aat',obj_divasl);


obj_divasl = obj_divasl.callCMF('white','all',1:11);
obj_divasl.dCMF.showResult('cbf',obj_divasl);
obj_divasl.dCMF.showResult('aat',obj_divasl);

% obj_divasl.dDT_CMF.showMRI('Div-TI, AAT',obj_divasl,'CMF');
%% display different slices and repetitions
% close all
slices =8;
rep = 1;
figure
for i = 1:11
%     current = obj_divasl.dDivasl(i).getMean('mri',:,:,slices,1);
    current = obj_divasl.dDivasl(i).getASL('mri',:,:,slices,rep);
    
    control = obj_divasl.dDivasl(i).dMRI(:,:,slices,(rep)*2);
    label = obj_divasl.dDivasl(i).dMRI(:,:,slices,(rep)*2+1);
    norm(current - (label-control))
    
    th = 5;
    als = @(control,th) (sqrt(sum(control(control>th).^2)/length(control(control>th))));
    energy1 = als(control(:),th);
    subplot(3,11,i),imshow(control,[]);title([num2str(energy1)]);
    energy2 = als(label(:),th);
    subplot(3,11,i+11),imshow(label,[]);title([num2str(energy2)]);
    energy3 = als(current(:),1);
    subplot(3,11,i+22),imshow(current,[-5,20]);title([num2str(energy3),',',num2str(energy3/((energy1+energy2)*2)*100)]);
end
return
% % Path to mask
% pathMask = 'C:\Users\lyu\Documents\MATLAB\myProject\pre-processing\gr_group01\su_subject01\structural';
% nameMask = 'c1rT13D_asl_file.nii';
% pathOpen = fullfile(pathMask,nameMask);
% dataMask = nii2mat(pathOpen);



% figure(1),imshow(obj_divasl.dDivasl(1).getASL('mri',:,:,8,10),[]);

% figure(2),imshow(obj_divasl.dDivasl(1).getMean('mri',:,:,8,1),[]);
%%
slices = 8;
for x = 5:64
    for y = 5:64
       plot( obj_divasl.getImg('divti',x,y,slices,1) );
       pause
    end
end

% figure(3),plot(obj_divasl.getImg('divti',[34,16,8,1]))
%%
